Fit density (also referred to as CPUE) model with environmental predictors and use that to calculate weighted mean dissolved oxygen, temperature and depth of Baltic cod
# Load libraries, install if needed
library(tidyverse); theme_set(theme_classic(base_size = 10))
library(readxl)
library(tidylog)
library(RCurl)
library(viridis)
library(RColorBrewer)
library(patchwork)
library(janitor)
library(icesDatras)
library(mapdata)
library(patchwork)
library(rgdal)
library(raster)
library(sf)
library(rgeos)
library(chron)
library(lattice)
library(ncdf4)
library(marmap)
library(rnaturalearth)
library(rnaturalearthdata)
library(mapplots)
library(EnvStats)
library(qwraps2)
library(bayesplot)
library(tidybayes)
library(brms)
#remotes::install_github("pbs-assess/sdmTMB")
library(sdmTMB) # To load entire cache in interactive r session, do: qwraps2::lazyload_cache_dir(path = "R/analysis/cod_flounder_diets_spatial_cache/html")
# Specify map ranges
ymin = 55; ymax = 58; xmin = 14; xmax = 20
map_data <- rnaturalearth::ne_countries(
scale = "medium",
returnclass = "sf", continent = "europe")
# Crop the polygon for plotting and efficiency:
# st_bbox(map_data) # find the rough coordinates
swe_coast <- suppressWarnings(suppressMessages(
st_crop(map_data,
c(xmin = xmin, ymin = ymin, xmax = xmax, ymax = ymax))))
# Transform our map into UTM 33 coordinates, which is the equal-area projection we fit in:
utm_zone33 <- 32633
swe_coast_proj <- sf::st_transform(swe_coast, crs = utm_zone33)
ggplot(swe_coast_proj) + geom_sf()
# Define plotting theme for main plot
theme_plot <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 8),
legend.text = element_text(size = 8),
legend.title = element_text(size = 8),
legend.position = "bottom",
legend.key.height = unit(0.2, "cm"),
legend.margin = margin(0, 0, 0, 0),
legend.box.margin = margin(-5, -5, -5, -5),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Define plotting theme for facet_wrap map with years
theme_facet_map <- function(base_size = 10, base_family = "") {
theme_light(base_size = 10, base_family = "") +
theme(
axis.text.x = element_text(angle = 90),
axis.text = element_text(size = 6),
strip.text = element_text(size = 8, colour = 'black', margin = margin()),
strip.background = element_rect(fill = "grey90")
)
}
# Make default base map plot
plot_map_raster <-
ggplot(swe_coast_proj) +
geom_sf(size = 0.3) +
labs(x = "Longitude", y = "Latitude") +
theme_facet_map(base_size = 14)
# Function to convert from lat lon to UTM
LongLatToUTM <- function(x, y, zone){
xy <- data.frame(ID = 1:length(x), X = x, Y = y)
coordinates(xy) <- c("X", "Y")
proj4string(xy) <- CRS("+proj=longlat +datum=WGS84") ## for example
res <- spTransform(xy, CRS(paste("+proj=utm +zone=",zone," ellps=WGS84",sep='')))
return(as.data.frame(res))
}
theme_clean <- function() {
theme_minimal(base_family = "Barlow Semi Condensed") +
theme(panel.grid.minor = element_blank(),
plot.title = element_text(family = "BarlowSemiCondensed-Bold"),
axis.title = element_text(family = "BarlowSemiCondensed-Medium"),
strip.text = element_text(family = "BarlowSemiCondensed-Bold",
size = rel(1), hjust = 0),
strip.background = element_rect(fill = "grey80", color = NA))
}
# These data are for total and prey specific stomach models
cod <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/cod_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> time_period = col_character(),
#> quarter_fact = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
cod <- cod %>%
mutate(year = as.integer(year),
quarter = as.factor(quarter),
depth2_sc = depth - mean(depth),
saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
filter(year > 2014) %>%
filter(!quarter == 2) %>%
drop_na(predfle_density_sc, predcod_density_sc) %>%
droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#> converted 'quarter' from double to factor (0 new NA)
#> new variable 'depth2_sc' (double) with 106 unique values and 0% NA
#> new variable 'saduria_entomon_per_mass' (double) with 1,709 unique values and 0% NA
#> new variable 'tot_prey_biom_per_mass' (double) with 7,202 unique values and 0% NA
#> new variable 'depth_sc' (double) with 106 unique values and 0% NA
#> filter: removed 4,968 rows (57%), 3,703 rows remaining
#> filter: removed 91 rows (2%), 3,612 rows remaining
#> drop_na: no rows removed
fle <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/fle_diet_analysis.csv") %>% dplyr::select(-X1)
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
fle <- fle %>%
mutate(year = as.integer(year),
quarter = as.factor(quarter),
depth2_sc = depth - mean(depth),
saduria_entomon_per_mass = saduria_entomon_tot/pred_weight_g,
tot_prey_biom_per_mass = tot_prey_biom/pred_weight_g,
depth_sc = (depth - mean(depth)) / sd(depth)) %>%
filter(!quarter == 2) %>%
drop_na(predfle_density_sc, predcod_density_sc) %>%
droplevels()
#> mutate: converted 'year' from double to integer (0 new NA)
#> converted 'quarter' from double to factor (0 new NA)
#> new variable 'depth2_sc' (double) with 81 unique values and 0% NA
#> new variable 'saduria_entomon_per_mass' (double) with 608 unique values and 0% NA
#> new variable 'tot_prey_biom_per_mass' (double) with 1,735 unique values and 0% NA
#> new variable 'depth_sc' (double) with 81 unique values and 0% NA
#> filter: removed 58 rows (3%), 2,190 rows remaining
#> drop_na: no rows removed
# These data are for the diet analysis (all prey groups)
cod_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/pca_cod_diet_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "COD") %>% filter(year > 2014) %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: removed 4,968 rows (57%), 3,703 rows remaining
#> filter: removed 91 rows (2%), 3,612 rows remaining
fle_prey <- readr::read_csv("/Users/maxlindmark/Desktop/R_STUDIO_PROJECTS/cod_interactions/data/pca_fle_diet_analysis.csv") %>% dplyr::select(-X1) %>% mutate(species = "FLE") %>% filter(!quarter == 2) %>% droplevels()
#> Warning: Missing column names filled in: 'X1' [1]
#>
#> ── Column specification ────────────────────────────────────────────────────────
#> cols(
#> .default = col_double(),
#> unique_pred_id = col_character(),
#> predator_spec = col_character(),
#> ices_rect = col_character(),
#> cruise = col_character()
#> )
#> ℹ Use `spec()` for the full column specifications.
#> mutate: new variable 'species' (character) with one unique value and 0% NA
#> filter: removed 58 rows (3%), 2,190 rows remaining
# And now read in pred_grid2 which has oxygen values at location and time and depth:
pred_grid2 <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/cod_condition/master/data/for_analysis/pred_grid2.csv") ### CHANGE TO OUR NEW PRED GRID
# Clean
pred_grid2 <- pred_grid2 %>%
mutate(year = as.integer(year)) %>%
drop_na(depth)
# Add ices_rect
pred_grid2$ices_rect <- ices.rect2(pred_grid2$lon, pred_grid2$lat)
# pred_grid2_q1 <- pred_grid2 %>% mutate(quarter = factor(1))
# pred_grid2_q4 <- pred_grid2 %>% mutate(quarter = factor(4))
# Plot data in space
plot_map_raster +
geom_point(data = fle_prey, aes(x = X * 1000, y = Y * 1000, color = "FLE"), alpha = 0.5) +
geom_point(data = cod_prey, aes(x = X * 1000, y = Y * 1000, color = "COD"), alpha = 0.5) +
theme_plot()
cod_prey <- cod_prey %>% filter(lat < 58)
#> filter: removed 19 rows (1%), 3,593 rows remaining
fle_prey <- fle_prey %>% filter(lat < 58)
#> filter: removed 10 rows (<1%), 2,180 rows remaining
# Reformat data to calculate Schoeners overlap index
# colnames(fle_prey)
fle_prey_long <- fle_prey %>%
pivot_longer(14:29) %>%
group_by(name, year, quarter, ices_rect) %>% # ices rect?
summarise(fle_diet = sum(value)) %>%
arrange(name, year, ices_rect) %>%
ungroup() %>%
group_by(year, quarter, ices_rect) %>%
mutate(fle_diet_tot = sum(fle_diet),
fle_diet_prop = fle_diet / fle_diet_tot) %>%
ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 2180x30, now 34880x16]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,296 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> mutate (grouped): new variable 'fle_diet_tot' (double) with 81 unique values and 0% NA
#> new variable 'fle_diet_prop' (double) with 445 unique values and 1% NA
#> ungroup: no grouping variables
fle_prey_wide <- fle_prey_long %>%
pivot_wider(names_from = name, values_from = fle_diet_prop, values_fill = 0)
#> pivot_wider: reorganized (name, fle_diet_prop) into (amphipoda_tot, bivalvia_tot, clupea_harengus_tot, clupeidae_tot, gadiformes_tot, …) [was 1296x7, now 526x21]
cod_prey_long <- cod_prey %>%
pivot_longer(14:29) %>%
group_by(name, year, quarter, ices_rect) %>%
summarise(cod_diet = sum(value)) %>%
arrange(name, year, ices_rect) %>%
ungroup() %>%
group_by(year, quarter, ices_rect) %>%
mutate(cod_diet_tot = sum(cod_diet),
cod_diet_prop = cod_diet / cod_diet_tot) %>%
ungroup()
#> pivot_longer: reorganized (amphipoda_tot, bivalvia_tot, clupeidae_tot, clupea_harengus_tot, gadus_morhua_tot, …) into (name, value) [was 3593x30, now 57488x16]
#> group_by: 4 grouping variables (name, year, quarter, ices_rect)
#> summarise: now 1,696 rows and 5 columns, 3 group variables remaining (name, year, quarter)
#> ungroup: no grouping variables
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> mutate (grouped): new variable 'cod_diet_tot' (double) with 106 unique values and 0% NA
#> new variable 'cod_diet_prop' (double) with 784 unique values and 0% NA
#> ungroup: no grouping variables
cod_prey_wide <- cod_prey_long %>%
pivot_wider(names_from = name, values_from = cod_diet_prop, values_fill = 0)
#> pivot_wider: reorganized (name, cod_diet_prop) into (amphipoda_tot, bivalvia_tot, clupea_harengus_tot, clupeidae_tot, gadiformes_tot, …) [was 1696x7, now 892x21]
# Calculate Schoener index
schoener <- left_join(cod_prey_long, fle_prey_long) %>%
drop_na(name) %>%
drop_na(fle_diet_prop) %>%
drop_na(cod_diet_prop) %>%
group_by(year, quarter, ices_rect) %>%
summarise(schoener = 1 - 0.5*(sum(abs(fle_diet_prop - cod_diet_prop)))) %>%
ungroup()
#> Joining, by = c("name", "year", "quarter", "ices_rect")
#> left_join: added 3 columns (fle_diet, fle_diet_tot, fle_diet_prop)
#> > rows only in x 432
#> > rows only in y ( 32)
#> > matched rows 1,264
#> > =======
#> > rows total 1,696
#> drop_na: no rows removed
#> drop_na: removed 448 rows (26%), 1,248 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 78 rows and 4 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables
# Calculate Levin niche index
levin <- left_join(cod_prey_long, fle_prey_long) %>%
drop_na(name) %>%
drop_na(fle_diet_prop) %>%
drop_na(cod_diet_prop) %>%
group_by(year, quarter, ices_rect) %>%
summarise(levin_cod = (1/(16-1)) * (((1/sum(cod_diet_prop^2))) - 1),
levin_fle = (1/(16-1)) * (((1/sum(fle_diet_prop^2))) - 1)) %>%
ungroup()
#> Joining, by = c("name", "year", "quarter", "ices_rect")
#> left_join: added 3 columns (fle_diet, fle_diet_tot, fle_diet_prop)
#> > rows only in x 432
#> > rows only in y ( 32)
#> > matched rows 1,264
#> > =======
#> > rows total 1,696
#> drop_na: no rows removed
#> drop_na: removed 448 rows (26%), 1,248 rows remaining
#> drop_na: no rows removed
#> group_by: 3 grouping variables (year, quarter, ices_rect)
#> summarise: now 78 rows and 5 columns, 2 group variables remaining (year, quarter)
#> ungroup: no grouping variables
# Merge the indicies
ind <- left_join(levin, schoener)
#> Joining, by = c("year", "quarter", "ices_rect")
#> left_join: added one column (schoener)
#> > rows only in x 0
#> > rows only in y ( 0)
#> > matched rows 78
#> > ====
#> > rows total 78
# Summarise cod and flounder data by ices_rect then add to diet data
colnames(cod)
#> [1] "FR_tot" "FR_sad"
#> [3] "FR_spr" "saduria_entomon_tot"
#> [5] "tot_prey_biom" "common_tot"
#> [7] "unique_pred_id" "year"
#> [9] "quarter" "time_period"
#> [11] "quarter_fact" "pred_weight_g"
#> [13] "pred_cm" "predator_spec"
#> [15] "predcod_density" "predfle_density"
#> [17] "predcod_density_sc" "predfle_density_sc"
#> [19] "depth" "X"
#> [21] "Y" "lat"
#> [23] "long" "ices_rect"
#> [25] "cruise" "depth2_sc"
#> [27] "saduria_entomon_per_mass" "tot_prey_biom_per_mass"
#> [29] "depth_sc"
cod$year_rect_id <- paste(cod$year, cod$quarter, cod$ices_rect, sep = "_")
dens_sum <- cod %>% group_by(year_rect_id) %>%
summarise(predfle_density_tot = sum(predfle_density),
predcod_density_tot = sum(predcod_density)) %>%
ungroup() %>%
mutate(predfle_density_tot_sc = (predfle_density_tot - mean(predfle_density_tot)) / sd(predfle_density_tot),
predcod_density_tot_sc = (predcod_density_tot - mean(predcod_density_tot)) / sd(predcod_density_tot))
#> group_by: one grouping variable (year_rect_id)
#> summarise: now 112 rows and 3 columns, ungrouped
#> ungroup: no grouping variables
#> mutate: new variable 'predfle_density_tot_sc' (double) with 110 unique values and 0% NA
#> new variable 'predcod_density_tot_sc' (double) with 110 unique values and 0% NA
ind$year_rect_id <- paste(ind$year, ind$quarter, ind$ices_rect, sep = "_")
ind <- left_join(ind, dens_sum)
#> Joining, by = "year_rect_id"
#> left_join: added 4 columns (predfle_density_tot, predcod_density_tot, predfle_density_tot_sc, predcod_density_tot_sc)
#> > rows only in x 0
#> > rows only in y (34)
#> > matched rows 78
#> > ====
#> > rows total 78
# Summarise depth from the prediction grid then add to diet data
pred_grid2_sum <- pred_grid2 %>%
group_by(ices_rect) %>%
summarise(mean_depth = mean(depth)) %>%
mutate(depth_sc = (mean_depth - mean(mean_depth)) / sd(mean_depth)) %>%
ungroup()
#> group_by: one grouping variable (ices_rect)
#> summarise: now 61 rows and 2 columns, ungrouped
#> mutate: new variable 'depth_sc' (double) with 61 unique values and 0% NA
#> ungroup: no grouping variables
ind <- left_join(ind, dplyr::select(pred_grid2_sum, ices_rect, depth_sc, mean_depth))
#> Joining, by = "ices_rect"
#> left_join: added 2 columns (depth_sc, mean_depth)
#> > rows only in x 0
#> > rows only in y (43)
#> > matched rows 78
#> > ====
#> > rows total 78
# Plot the indicies
ggplot(ind) +
geom_jitter(aes(factor(year), schoener),
alpha = 0.8, width = 0.2, height = 0, size = 2)
ggplot(ind) +
geom_jitter(aes(factor(year), levin_cod, color = "cod"),
alpha = 0.8, width = 0.2, height = 0, size = 2) +
geom_jitter(aes(factor(year), levin_fle, color = "fle"),
alpha = 0.8, width = 0.2, height = 0, size = 2) +
scale_color_brewer(palette = "Dark2")
# Against response variables
# Schoener
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = schoener)) +
geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
ggplot(ind, aes(x = predfle_density_tot, y = schoener)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
ggplot(ind, aes(x = predcod_density_tot, y = schoener)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
# Now Levin
# Cod
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = levin_cod)) +
geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
ggplot(ind, aes(x = predfle_density_tot, y = levin_cod)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
ggplot(ind, aes(x = predcod_density_tot, y = levin_cod)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
# Flounder
ggplot(ind, aes(x = predfle_density_tot + predcod_density_tot, y = levin_fle)) +
geom_point() + stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
ggplot(ind, aes(x = predfle_density_tot, y = levin_fle)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
ggplot(ind, aes(x = predcod_density_tot, y = levin_fle)) + geom_point() +
stat_smooth(method = "gam", formula = y ~ s(x, k = 3)) + facet_wrap(~quarter)
# Add refs here or in text already! THis is complementary to prey weights in stomach
# Do cumsum plots
# Now split by size-classes, check Orior co author paper for which sizes to split
# Think about which model to fit to these indices?
# ---- brms with random rectangle and fixed year? Evaluate differences in intercept and plot covariates?
# ... Then when I have size-group? I get like 5 "species". Maybe do pca on those?
# Do that after resolving old comments, adding new diet data and correctly adding quarter and new cpue data!!!
# Quickly check data to determine which distribution to use
ind %>%
ungroup() %>%
count(schoener == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#> `schoener == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 75 0.962
#> 2 TRUE 3 0.0385
ind %>%
ungroup() %>%
count(levin_cod == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#> `levin_cod == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 77 0.987
#> 2 TRUE 1 0.0128
ind %>%
ungroup() %>%
count(levin_fle == 0) %>%
mutate(prop = n / sum(n))
#> ungroup: no grouping variables
#> count: now 2 rows and 2 columns, ungrouped
#> mutate: new variable 'prop' (double) with 2 unique values and 0% NA
#> # A tibble: 2 × 3
#> `levin_fle == 0` n prop
#> <lgl> <int> <dbl>
#> 1 FALSE 75 0.962
#> 2 TRUE 3 0.0385
# Fit beta models, so few zeroes!
ind %>% arrange(schoener) %>% dplyr::select(schoener)
#> # A tibble: 78 × 1
#> schoener
#> <dbl>
#> 1 0
#> 2 0
#> 3 0
#> 4 0.000490
#> 5 0.00127
#> 6 0.00134
#> 7 0.00149
#> 8 0.00209
#> 9 0.00240
#> 10 0.00261
#> # … with 68 more rows
ind %>% arrange(levin_cod) %>% dplyr::select(levin_cod)
#> # A tibble: 78 × 1
#> levin_cod
#> <dbl>
#> 1 0
#> 2 0.00173
#> 3 0.00265
#> 4 0.00493
#> 5 0.00801
#> 6 0.0132
#> 7 0.0163
#> 8 0.0179
#> 9 0.0191
#> 10 0.0270
#> # … with 68 more rows
ind %>% arrange(levin_fle) %>% dplyr::select(levin_fle)
#> # A tibble: 78 × 1
#> levin_fle
#> <dbl>
#> 1 0
#> 2 0
#> 3 0
#> 4 0.000321
#> 5 0.000737
#> 6 0.00191
#> 7 0.00206
#> 8 0.00304
#> 9 0.00779
#> 10 0.00811
#> # … with 68 more rows
ind <- ind %>% mutate(schoener2 = ifelse(schoener == 0, 0.0001, schoener),
levin_cod2 = ifelse(levin_cod == 0, 0.0001, levin_cod),
levin_fle2 = ifelse(levin_fle == 0, 0.0001, levin_fle),
year_f = as.factor(year),
quarter_f = as.factor(quarter),
ices_rect_f = as.factor(ices_rect))
#> mutate: new variable 'schoener2' (double) with 76 unique values and 0% NA
#> new variable 'levin_cod2' (double) with 78 unique values and 0% NA
#> new variable 'levin_fle2' (double) with 76 unique values and 0% NA
#> new variable 'year_f' (factor) with 6 unique values and 0% NA
#> new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> new variable 'ices_rect_f' (factor) with 18 unique values and 0% NA
brms models with densities as covariates to diversity and overlap indices# All covariates
m_schoen_beta_full <- brm(
bf(schoener2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f + quarter_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Trying to compile a simple C file
#> Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
#> clang -mmacosx-version-min=10.13 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Users/maxlindmark/Library/R/4.0/library/Rcpp/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/unsupported" -I"/Users/maxlindmark/Library/R/4.0/library/BH/include" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/src/" -I"/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/" -I"/Users/maxlindmark/Library/R/4.0/library/RcppParallel/include/" -I"/Users/maxlindmark/Library/R/4.0/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DBOOST_NO_AUTO_PTR -include '/Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/usr/local/include -fPIC -Wall -g -O2 -c foo.c -o foo.o
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:88:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:1: error: unknown type name 'namespace'
#> namespace Eigen {
#> ^
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:628:16: error: expected ';' after top level declarator
#> namespace Eigen {
#> ^
#> ;
#> In file included from <built-in>:1:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/StanHeaders/include/stan/math/prim/mat/fun/Eigen.hpp:13:
#> In file included from /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Dense:1:
#> /Users/maxlindmark/Library/R/4.0/library/RcppEigen/include/Eigen/Core:96:10: fatal error: 'complex' file not found
#> #include <complex>
#> ^~~~~~~~~
#> 3 errors generated.
#> make: *** [foo.o] Error 1
#> Start sampling
plot(m_schoen_beta_full)
conditional_effects(m_schoen_beta_full)
loo_m_schoen_beta_full <- loo(m_schoen_beta_full, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta_full)
# Only density covariates
m_schoen_beta_dens <- brm(
bf(schoener2 ~ 0 + year_f + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Start sampling
plot(m_schoen_beta_dens)
conditional_effects(m_schoen_beta_dens)
loo_m_schoen_beta_dens <- loo(m_schoen_beta_dens, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta_dens)
# Simplest model
m_schoen_beta <- brm(
bf(schoener2 ~ 0 + year_f + (1|ices_rect_f),
phi ~ 0 + year_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 4, iter = 4000, warmup = 1000, cores = 4, control = list(adapt_delta = 0.95))
#> Compiling Stan program...
#> Start sampling
plot(m_schoen_beta)
conditional_effects(m_schoen_beta)
loo_m_schoen_beta <- loo(m_schoen_beta, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(loo_m_schoen_beta)
loo_compare(loo_m_schoen_beta, loo_m_schoen_beta_dens, loo_m_schoen_beta_full)
#> elpd_diff se_diff
#> m_schoen_beta 0.0 0.0
#> m_schoen_beta_dens -1.0 1.4
#> m_schoen_beta_full -4.3 2.6
# Plot models (working with the m_schoen_beta_dens)
summary(m_schoen_beta_dens)
#> Family: beta
#> Links: mu = logit; phi = log
#> Formula: schoener2 ~ 0 + year_f + predfle_density_tot_sc + predcod_density_tot_sc + (1 | ices_rect_f)
#> phi ~ 0 + year_f
#> Data: ind (Number of observations: 78)
#> Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
#> total post-warmup draws = 12000
#>
#> Group-Level Effects:
#> ~ices_rect_f (Number of levels: 18)
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept) 0.19 0.16 0.01 0.58 1.00 4846 5380
#>
#> Population-Level Effects:
#> Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> year_f2015 -2.74 0.46 -3.59 -1.75 1.00 10133
#> year_f2016 -2.91 0.30 -3.48 -2.29 1.00 11525
#> year_f2017 -2.57 0.36 -3.24 -1.83 1.00 11899
#> year_f2018 -2.14 0.51 -3.08 -1.08 1.00 10237
#> year_f2019 -1.58 0.51 -2.55 -0.55 1.00 12042
#> year_f2020 -1.46 0.42 -2.26 -0.60 1.00 11745
#> predfle_density_tot_sc -0.25 0.16 -0.57 0.06 1.00 10524
#> predcod_density_tot_sc 0.18 0.15 -0.11 0.48 1.00 10331
#> phi_year_f2015 2.36 0.64 0.95 3.46 1.00 9842
#> phi_year_f2016 2.30 0.38 1.49 2.97 1.00 11640
#> phi_year_f2017 1.82 0.42 0.91 2.56 1.00 12310
#> phi_year_f2018 1.21 0.54 0.02 2.17 1.00 10486
#> phi_year_f2019 0.67 0.51 -0.40 1.58 1.00 11405
#> phi_year_f2020 1.46 0.50 0.40 2.35 1.00 12856
#> Tail_ESS
#> year_f2015 7731
#> year_f2016 7531
#> year_f2017 8457
#> year_f2018 7513
#> year_f2019 8082
#> year_f2020 8680
#> predfle_density_tot_sc 9258
#> predcod_density_tot_sc 8973
#> phi_year_f2015 7032
#> phi_year_f2016 7508
#> phi_year_f2017 7677
#> phi_year_f2018 7340
#> phi_year_f2019 8729
#> phi_year_f2020 8334
#>
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
# Working with the posterior
posterior_beta <- m_schoen_beta_dens %>%
gather_draws(`b_.*`, regex = TRUE) %>%
mutate(component = ifelse(str_detect(.variable, "phi_"), "Precision", "Mean"),
intercept = str_detect(.variable, "Intercept"))
ggplot(posterior_beta, aes(x = .value, y = fct_rev(.variable), fill = component)) +
geom_vline(xintercept = 0) +
stat_halfeye(aes(slab_alpha = intercept), alpha = 0.5,
.width = c(0.8, 0.95), point_interval = "median_hdi") +
scale_fill_brewer(palette = "Dark2") +
scale_slab_alpha_discrete(range = c(1, 0.4)) +
guides(fill = "none", slab_alpha = "none") +
labs(x = "Coefficient", y = "Variable") +
facet_wrap(vars(component), ncol = 1, scales = "free_y") +
theme_plot() +
NULL
# Selected parameters
# https://www.sciencedirect.com/topics/mathematics/logit-scale
posterior_beta %>%
filter(component == "Mean" & .variable %in% c("b_predcod_density_tot_sc", "b_predfle_density_tot_sc")) %>%
ggplot(., aes(x = plogis(.value), y = .variable, fill = .variable)) +
geom_vline(xintercept = 0) +
stat_halfeye(alpha = 0.5, .width = c(0.8, 0.95)) +
scale_fill_brewer(palette = "Dark2") +
labs(x = "Coefficient", y = "Variable") +
theme_plot() +
guides(fill = FALSE) +
NULL
#> filter (grouped): removed 144,000 rows (86%), 24,000 rows remaining
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
# Prediction! (non dens model)
beta_bayes_pred_yr <- m_schoen_beta %>%
epred_draws(newdata = tibble(year_f = as.factor(c(2015:2020))),
re_formula = NA)
ggplot(beta_bayes_pred_yr, aes(x = .epred, y = year_f, color = year_f, fill = year_f)) +
stat_halfeye(.width = c(0.8, 0.95), alpha = 0.7) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
guides(fill = FALSE, color = FALSE) +
coord_flip(xlim = c(0, 0.4)) +
labs(x = "Predicted Schoener's overlap index", y = NULL) +
theme_plot()
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
beta_bayes_pred_re <- m_schoen_beta %>%
epred_draws(newdata = tibble(expand.grid(year_f = as.factor(c(2015:2020)),
ices_rect_f = as.factor(unique(ind$ices_rect)))))
ggplot(beta_bayes_pred_re, aes(x = .epred, y = ices_rect_f, color = ices_rect_f, fill = ices_rect_f)) +
stat_halfeye(.width = c(0.8, 0.95), alpha = 0.7) +
scale_fill_viridis(discrete = TRUE) +
scale_color_viridis(discrete = TRUE) +
guides(fill = FALSE, color = FALSE) +
coord_cartesian(xlim = c(0, 0.4)) +
#coord_flip(xlim = c(0, 0.4)) +
labs(x = "Predicted Schoener's overlap index", y = NULL) +
theme_plot()
#> Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
#> "none")` instead.
# Cod Levin index: all covariates
m_levin_beta_full_cod <- brm(
bf(levin_cod2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f + quarter_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 2, iter = 2000, warmup = 1000, cores = 2, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
plot(m_levin_beta_full_cod)
conditional_effects(m_levin_beta_full_cod)
loo_m_schoen_beta_full <- loo(m_levin_beta_full_cod, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(m_levin_beta_full_cod)
m_levin_beta_full_fle <- brm(
bf(levin_fle2 ~ 0 + year_f + quarter_f + depth_sc + predfle_density_tot_sc + predcod_density_tot_sc + (1|ices_rect_f),
phi ~ 0 + year_f + quarter_f),
data = ind, family = brms::Beta(link = "logit", link_phi = "log"), save_pars = save_pars(all = TRUE),
#backend = "cmdstanr",
chains = 2, iter = 2000, warmup = 1000, cores = 2, control = list(adapt_delta = 0.99))
#> Compiling Stan program...
#> Start sampling
plot(m_levin_beta_full_fle)
conditional_effects(m_levin_beta_full_fle)
loo_m_schoen_beta_full <- loo(m_levin_beta_full_fle, moment_match = TRUE)
#> Warning: Some Pareto k diagnostic values are slightly high. See help('pareto-k-diagnostic') for details.
plot(m_levin_beta_full_fle)
sdmTMB models with densities as covariates to stomach content data. First make mesh# Cod
pcod_spde <- make_mesh(cod, c("X", "Y"), n_knots = 150, type = "kmeans", seed = 42)
plot(pcod_spde)
pfle_spde <- make_mesh(fle, c("X", "Y"), n_knots = 70, type = "kmeans", seed = 42)
plot(pfle_spde)
# Model total prey biomass
# Cod
mcod <- sdmTMB(
data = cod,
formula = tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc + s(depth_sc),
time = "year", fields = "IID", spatial_only = FALSE, spde = pcod_spde, family = tweedie())
#> Registered S3 methods overwritten by 'lme4':
#> method from
#> cooks.distance.influence.merMod car
#> influence.merMod car
#> dfbeta.influence.merMod car
#> dfbetas.influence.merMod car
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k` since the degree of wiggliness
#> is determined by the data.
tidy(mcod)
#> term estimate std.error
#> 1 quarter1 -4.81710803 0.18724912
#> 2 quarter4 -4.66825835 0.16127572
#> 3 as.factor(year)2016 0.24009064 0.19659901
#> 4 as.factor(year)2017 0.36443934 0.19436610
#> 5 as.factor(year)2018 0.42481751 0.21019866
#> 6 as.factor(year)2019 0.23688896 0.23211806
#> 7 as.factor(year)2020 0.47897776 0.23197072
#> 8 predfle_density_sc -0.02663818 0.02644890
#> 9 predcod_density_sc 0.09099320 0.03395342
summary(mcod)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc +
#> Formula: predcod_density_sc + s(depth_sc)
#> Time column: "year"
#> SPDE: pcod_spde
#> Data: cod
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> quarter1 -4.82 0.19
#> quarter4 -4.67 0.16
#> as.factor(year)2016 0.24 0.20
#> as.factor(year)2017 0.36 0.19
#> as.factor(year)2018 0.42 0.21
#> as.factor(year)2019 0.24 0.23
#> as.factor(year)2020 0.48 0.23
#> predfle_density_sc -0.03 0.03
#> predcod_density_sc 0.09 0.03
#>
#> Dispersion parameter: 0.52
#> Tweedie p: 1.74
#> Matern range: 9.74
#> Spatial SD: 0.23
#> Spatiotemporal SD: 0.68
#> ML criterion at convergence: -9961.622
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
cod$resids_mcod <- residuals(mcod) # randomized quantile residuals
qqnorm(cod$resids_mcod); abline(a = 0, b = 1)
# Flounder
mfle <- sdmTMB(
data = fle,
formula = tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc + s(depth_sc),
time = "year", fields = "IID", spatial_only = FALSE, spde = pfle_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k` since the degree of wiggliness
#> is determined by the data.
tidy(mfle)
#> term estimate std.error
#> 1 quarter1 -5.038155329 0.22043580
#> 2 quarter4 -4.544660978 0.18833024
#> 3 as.factor(year)2016 -0.020008080 0.24180860
#> 4 as.factor(year)2017 0.378980202 0.23768869
#> 5 as.factor(year)2018 0.593044431 0.28817619
#> 6 as.factor(year)2019 0.584962410 0.25846307
#> 7 as.factor(year)2020 0.995539716 0.25328069
#> 8 predfle_density_sc -0.043565572 0.03287689
#> 9 predcod_density_sc 0.002175145 0.03343344
summary(mfle)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: tot_prey_biom_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc +
#> Formula: predcod_density_sc + s(depth_sc)
#> Time column: "year"
#> SPDE: pfle_spde
#> Data: fle
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> quarter1 -5.04 0.22
#> quarter4 -4.54 0.19
#> as.factor(year)2016 -0.02 0.24
#> as.factor(year)2017 0.38 0.24
#> as.factor(year)2018 0.59 0.29
#> as.factor(year)2019 0.58 0.26
#> as.factor(year)2020 1.00 0.25
#> predfle_density_sc -0.04 0.03
#> predcod_density_sc 0.00 0.03
#>
#> Dispersion parameter: 0.12
#> Tweedie p: 1.50
#> Matern range: 29.17
#> Spatial SD: 0.30
#> Spatiotemporal SD: 0.54
#> ML criterion at convergence: -4721.832
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
fle$resids_mfle <- residuals(mfle) # randomized quantile residuals
qqnorm(fle$resids_mfle); abline(a = 0, b = 1)
# Model Saduria contents
# Cod
mcodsad <- sdmTMB(
data = cod,
formula = saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc + s(depth_sc),
time = "year", fields = "IID", spatial_only = FALSE, spde = pcod_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k` since the degree of wiggliness
#> is determined by the data.
tidy(mcodsad)
#> term estimate std.error
#> 1 quarter1 -9.02384114 0.7872885
#> 2 quarter4 -9.02682697 0.7395658
#> 3 as.factor(year)2016 -1.07056118 0.7875507
#> 4 as.factor(year)2017 -1.45889570 0.7842763
#> 5 as.factor(year)2018 -0.21692129 0.8426281
#> 6 as.factor(year)2019 -1.04243883 0.9148954
#> 7 as.factor(year)2020 0.06204430 0.8954815
#> 8 predfle_density_sc -0.10384084 0.1238834
#> 9 predcod_density_sc -0.04374943 0.1062921
summary(mcodsad)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc +
#> Formula: predcod_density_sc + s(depth_sc)
#> Time column: "year"
#> SPDE: pcod_spde
#> Data: cod
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> quarter1 -9.02 0.79
#> quarter4 -9.03 0.74
#> as.factor(year)2016 -1.07 0.79
#> as.factor(year)2017 -1.46 0.78
#> as.factor(year)2018 -0.22 0.84
#> as.factor(year)2019 -1.04 0.91
#> as.factor(year)2020 0.06 0.90
#> predfle_density_sc -0.10 0.12
#> predcod_density_sc -0.04 0.11
#>
#> Dispersion parameter: 0.17
#> Tweedie p: 1.49
#> Matern range: 41.79
#> Spatial SD: 2.00
#> Spatiotemporal SD: 1.53
#> ML criterion at convergence: -1022.716
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
cod$resids_mcodsad <- residuals(mcodsad) # randomized quantile residuals
qqnorm(cod$resids_mcodsad); abline(a = 0, b = 1)
# Flounder
mflesad <- sdmTMB(
data = fle,
formula = saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc + predcod_density_sc + s(depth_sc),
time = "year", fields = "IID", spatial_only = FALSE, spde = pfle_spde, family = tweedie())
#> Warning: Detected a `s()` smoother. Smoothers are penalized in sdmTMB as
#> of version 0.0.19, but used to be unpenalized.
#> You no longer need to specify `k` since the degree of wiggliness
#> is determined by the data.
tidy(mflesad)
#> term estimate std.error
#> 1 quarter1 -8.99399061 0.98263924
#> 2 quarter4 -7.87852100 0.94515317
#> 3 as.factor(year)2016 1.00546603 0.80258800
#> 4 as.factor(year)2017 0.99471876 0.79567385
#> 5 as.factor(year)2018 1.89586027 0.91527579
#> 6 as.factor(year)2019 0.52603927 0.85483508
#> 7 as.factor(year)2020 2.15101984 0.82404480
#> 8 predfle_density_sc 0.02410729 0.07481394
#> 9 predcod_density_sc -0.17705506 0.07777392
summary(mflesad)
#> Spatial model fit by ML ['sdmTMB']
#> Formula: saduria_entomon_per_mass ~ 0 + quarter + as.factor(year) + predfle_density_sc +
#> Formula: predcod_density_sc + s(depth_sc)
#> Time column: "year"
#> SPDE: pfle_spde
#> Data: fle
#> Family: tweedie(link = 'log')
#> coef.est coef.se
#> quarter1 -8.99 0.98
#> quarter4 -7.88 0.95
#> as.factor(year)2016 1.01 0.80
#> as.factor(year)2017 0.99 0.80
#> as.factor(year)2018 1.90 0.92
#> as.factor(year)2019 0.53 0.85
#> as.factor(year)2020 2.15 0.82
#> predfle_density_sc 0.02 0.07
#> predcod_density_sc -0.18 0.08
#>
#> Dispersion parameter: 0.14
#> Tweedie p: 1.47
#> Matern range: 77.53
#> Spatial SD: 2.11
#> Spatiotemporal SD: 1.18
#> ML criterion at convergence: -1467.275
#>
#> See ?tidy.sdmTMB to extract these values as a data frame.
fle$resids_mflesad <- residuals(mflesad) # randomized quantile residuals
qqnorm(fle$resids_mflesad); abline(a = 0, b = 1)
# Flounder density
nd_cod_fle <- data.frame(predfle_density_sc = seq(min(cod$predfle_density_sc),
max(cod$predfle_density_sc), length.out = 100))
nd_cod_fle$year <- 2018L
nd_cod_fle$predcod_density_sc <- 0
nd_cod_fle$depth_sc <- 0
nd_cod_fle$quarter <- factor(4)
sad_pred_cod_fle <- predict(mcodsad, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)
tot_pred_cod_fle <- predict(mcod, newdata = nd_cod_fle, se_fit = TRUE, re_form = NA)
# Cod density
nd_cod_cod <- data.frame(predcod_density_sc = seq(min(cod$predcod_density_sc),
max(cod$predcod_density_sc), length.out = 100))
nd_cod_cod$year <- 2018L
nd_cod_cod$predfle_density_sc <- 0
nd_cod_cod$depth_sc <- 0
nd_cod_cod$quarter <- factor(4)
sad_pred_cod_cod <- predict(mcodsad, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)
tot_pred_cod_cod <- predict(mcod, newdata = nd_cod_cod, se_fit = TRUE, re_form = NA)
# Flounder density
nd_fle_fle <- data.frame(predfle_density_sc = seq(min(fle$predfle_density_sc),
max(fle$predfle_density_sc), length.out = 100))
nd_fle_fle$year <- 2018L
nd_fle_fle$predcod_density_sc <- 0
nd_fle_fle$depth_sc <- 0
nd_fle_fle$quarter <- factor(4)
sad_pred_fle_fle <- predict(mflesad, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)
tot_pred_fle_fle <- predict(mfle, newdata = nd_fle_fle, se_fit = TRUE, re_form = NA)
# Cod density
nd_fle_cod <- data.frame(predcod_density_sc = seq(min(fle$predcod_density_sc),
max(fle$predcod_density_sc), length.out = 100))
nd_fle_cod$year <- 2018L
nd_fle_cod$predfle_density_sc <- 0
nd_fle_cod$depth_sc <- 0
nd_fle_cod$quarter <- factor(4)
sad_pred_fle_cod <- predict(mflesad, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)
tot_pred_fle_cod <- predict(mfle, newdata = nd_fle_cod, se_fit = TRUE, re_form = NA)
# Saduria models
p1 <- ggplot(sad_pred_fle_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
#coord_cartesian(ylim = c(0, 0.015), expand = 0) +
labs(y = "Saduria in flounder stomach [g/g]", x = "")
p2 <- ggplot(sad_pred_fle_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(0, 0.015), expand = 0) +
labs(y = "", x = "")
p3 <- ggplot(sad_pred_cod_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "Saduria in cod stomach [g/g]", x = "Scaled flounder density")
p4 <- ggplot(sad_pred_cod_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
labs(y = "", x = "Scaled cod density")
p1 + p2 + p3 + p4 + plot_layout(ncol = 2)
ggsave("figures/marginal_effects_saduria.png", width = 6.5, height = 6.5, dpi = 600)
# Total stomach content models
p5 <- ggplot(tot_pred_fle_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(0, 0.033), expand = 0) +
labs(y = "Prey biomass in flounder stomach [g/g]", x = "")
p6 <- ggplot(tot_pred_fle_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(0, 0.033), expand = 0) +
labs(y = "", x = "")
p7 <- ggplot(tot_pred_cod_fle, aes(predfle_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(0, 0.028), expand = 0) +
labs(y = "Prey biomass in cod stomach [g/g]", x = "Scaled flounder density")
p8 <- ggplot(tot_pred_cod_cod, aes(predcod_density_sc, exp(est),
ymin = exp(est - 1.96 * est_se), ymax = exp(est + 1.96 * est_se))) +
geom_line() + geom_ribbon(alpha = 0.4) +
coord_cartesian(ylim = c(0, 0.028), expand = 0) +
labs(y = "", x = "Scaled cod density")
p5 + p6 + p7 + p8 + plot_layout(ncol = 2)
ggsave("figures/marginal_effects_total.png", width = 6.5, height = 6.5, dpi = 600)
knitr::knit_exit()